This notebook replicates key aspects of Young, Old, Conservative, and Bold by Gârleanu and Panageas published in the Journal of Political Economy in 2015.
The road map for this notebook is:
The model evolves around two types of households, both have Epstein-Zin-Weil preferences. $A$ is bold, takes on more risk and reacts more to a change in the interest rate. On the contrary, $B$ is conservative. Any household has an income that is hump-shaped in its age and death occurs according to a Poisson process with constant intensity, $\pi$, and they have no bequest motives.
Production, $Y_t$, follows a Brownian motion with drift and takes place in a firm. A constant fraction, $\omega$, of production is paid as income to the households, and the rest is distributed as dividends. There is in total one share of this firm and its price is $S_t$. There is also a risk-free bond in zero net supply that has a return $r_t$.
The equilibrium is described by the process for the consumption share of type-A households, $X_t$, which is the only state variable:
$$dX_t = \mu_X(X_t) dt + \sigma_X(X_t)dB_t$$
The above drift and volatility are functions of the solution objects, $g^A$, $g^B$, $\phi^1$, and $\phi^2$. $g^i$ is the consumption-to-wealth ratio and $\phi=\phi^1+\phi^2$ is the ratio of the present value of earnings at birth to aggregate consumption. There are two $\phi^j$ functions because the hump-shape of life-cycle income is parametrized by the sum of two exponential functions.
The equilibrium is described by the following system of coupled, second-order ODEs:
\begin{equation} 0=\frac{\sigma_X^2}{2}\frac{d^2\phi^j}{dX_t^2}+\frac{d\phi^j}{dX_t}\big(\mu_X + \sigma_X(\sigma_Y-\kappa)\big) + \phi^j ( \mu_Y - r - \pi - \delta_j -\sigma_Y \kappa) + B_j \omega \end{equation}
\begin{equation} 0=\frac{\sigma_X^2}{2} M_1^i \Biggm ((M_1^i-1) \biggm(\frac{\frac{dg^i}{dX_t}}{g^i} \biggm)^2+\frac{\frac{d^2g^i}{dX_t^2}}{g^i} \Biggm ) + M_1^i\frac{\frac{dg^i}{dX_t}}{g^i}\big(\mu_X-M_2^i\sigma_X\kappa\big) + \\ \Biggm ( \frac{\kappa^2(X_t)}{2}M_2^i(M_2^i-1)-M_2^i\big(r(X_t)+\pi\big)-M_1^ig^i+\frac{\Xi_3^i}{\gamma^i}\Biggm) \end{equation}
where $\gamma^i$ is relative risk aversion, $\alpha^i = 1 - \frac{1}{\psi^i}$, $\psi^i$ is the elasticity of intertemporal substitution, $\kappa$ is the Sharpe ratio, $\delta^j$ and $B_j$ parametrize the life-cycle income, and $\omega$ is the fraction of production that is paid as income. All other objects are constants that are defined below.
Mind: Equation (A.22) in the paper, which is equation (3) here, has a typo at the second derivative of $g^i$.
The boundary conditions are, for $X_t = 0$ or $X_t = 1$: $$0 = \frac{d\phi^j}{dX_t} \mu_X + \phi^j\big(\mu_Y - r - \pi - \delta_j - \sigma_Y\kappa\big) + B_j \omega$$
$$0 = M_1^i \frac{1}{g^i}\frac{dg^i}{dX_t}\mu_X + \Big(\frac{\kappa(X_t)^2}{2}M_2^i(M_2^i - 1) - M_2^i(r(X_t) + \pi) - M_1^i g^i + \frac{\Xi_3^i}{\gamma^i}\Big)$$
It turns out after a lot of sleuthing, that formulating the system of equations with $g^i$ is much less stable than defining $p^i = \frac{1}{g^i}$ and translating it. The fractions in equation (2) above and the definitions below are substituted with:
$$\frac{\frac{dg^i}{dX_t}}{g^i} = -\frac{\frac{dp^i}{dX_t}}{p^i}$$ and $$\frac{\frac{d^2g^i}{dX_t^2}}{g^i} = 2\Bigg(\frac{\frac{dp^i}{dX_t}}{p^i}\Bigg)^2-\frac{\frac{d^2p^i}{dX_t^2}}{p^i}$$
The auxiliary functions that appear in the ODEs are defined the following way:
$$\sigma_{X}(X_{t})=\frac{X_{t}\big(\Gamma(X_{t})-\gamma^{A}\big)}{\frac{\Gamma(X_{t})}{\gamma^{B}}X_{t}(1-X_{t})\Big[\frac{1-\gamma^{A}-\alpha^{A}}{\alpha^{A}}\frac{g^{A'}}{g^{A}}-\frac{1-\gamma^{B}-\alpha^{B}}{\alpha^{B}}\frac{g^{B'}}{g^{B}}\Big]+\gamma^{A}}\sigma_{Y}$$
σ_X(X_t, pA, pAprime, pB, pBprime) = σ_Y * X_t * (Γ(X_t) - γA) /
(Γ(X_t)/γB * X_t * (1 - X_t) *
((1 - γA - αA) / αA * (-pAprime / pA) - (1 - γB - αB) / αB * (- pBprime / pB))
+ γA);
$$\mu_X(X_t) = X_t\Big[\frac{r(X_t)-\rho}{1-\alpha^A} + n^A(X_t)-\pi-\mu_Y\Big] + v^A\pi\beta^A(X_t)-\sigma_Y \sigma_X(X_t)$$
μ_X(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) = X_t *
((r(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) - ρ) / (1 - αA) + nA(X_t, pA, pAprime, pB, pBprime) - π - μ_Y) +
vA * π * β(pA, ϕ1, ϕ2) - σ_Y * σ_X(X_t, pA, pAprime, pB, pBprime);
$$\kappa(X_t) = \Gamma(X_t)\sigma_Y + \sum_i\omega^i(X_t)\Big(\frac{1-\gamma^i-\alpha^i}{\alpha^i}\Big) \frac{g^{i'}}{g^i}\sigma_X(X_t)$$
κ(X_t, pA, pAprime, pB, pBprime) = Γ(X_t) * σ_Y +
ωA(X_t) * (1-γA-αA) / αA * (-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime) +
ωB(X_t) * (1-γB-αB) / αB * (-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime);
$$r(X_{t})=\rho+\frac{1}{\Theta(X_{t})}\left\{ \mu_{Y}-\pi\Big(\sum_{i}v^{i}\beta^{i}(X_{t})-1\Big)\right\} -\frac{1}{\Theta(X_{t})}\sum_{i}X_{t}^{i}n^{i}(X_{t})$$
r(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) = ρ +
1 / Θ(X_t) * (μ_Y - π * (vA * β(pA, ϕ1, ϕ2) + vB * β(pB, ϕ1, ϕ2) - 1)) -
1 / Θ(X_t) * (X_t * nA(X_t, pA, pAprime, pB, pBprime) + (1 - X_t) * nB(X_t, pA, pAprime, pB, pBprime));
Mind: The second line in equation (A.16) is not correct, see this issue on GitHub!
Properly derived the definition of nA and nB is:
$$n^{i}(X_{t})=\frac{2-\alpha^{i}}{2\gamma^{i}(1-\alpha^{i})}\kappa^{2}(X_{t})+\frac{\alpha^{i}+\gamma^{i}-1}{2\gamma^{i}\alpha^{i}}\Bigg(\frac{g^{i'}}{g^{i}}\sigma_{X}(X_{t})\Bigg)^{2}+\frac{\alpha^{i}+\gamma^{i}-1}{\alpha^{i}\gamma^{i}}\Big(\frac{g^{i'}}{g^{i}}\sigma_{X}(X_{t})\Big)\kappa(X_{t}) $$
nA(X_t, pA, pAprime, pB, pBprime) = (2 - αA) / (2 * γA * (1-αA)) * κ(X_t, pA, pAprime, pB, pBprime)^2 +
(αA + γA - 1) / (2 * γA * αA) * ((-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime))^2 +
(αA + γA - 1) / (γA * αA) * ((-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime)) * κ(X_t, pA, pAprime, pB, pBprime)
nB(X_t, pA, pAprime, pB, pBprime) = (2 - αB) / (2 * γB * (1-αB)) * κ(X_t, pA, pAprime, pB, pBprime)^2 +
(αB + γB - 1) / (2 * γB * αB) * ((-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime))^2 +
(αB + γB - 1) / (γB * αB) * ((-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime)) * κ(X_t, pA, pAprime, pB, pBprime)
$$\beta^i(X_t)=g^i(X_t)\underset{\phi^2(X_t) + \phi^2(X_t)}{\underbrace{\phi(X_t)}}$$
β(p, ϕ1, ϕ2) = (1 / p) * (ϕ1 + ϕ2);
$$X_t^A = X_t\quad\text{and}\quad X_t^B = 1-X_t$$
$$\Gamma(X_t) = \frac{1}{\sum_i \frac{X_t^i}{\gamma^i}}$$
Γ(X_t) = 1 / (X_t / γA + (1 - X_t) / γB);
$$\Theta(X_t) = \sum_i \frac{X_t^i}{1-\alpha^i}$$
Θ(X_t) = X_t / (1 - αA) + (1 - X_t) / (1 - αB);
$$\omega^i(X_t) = X_t^i \frac{\Gamma(X_t)}{\gamma^i}$$
ωA(X_t) = X_t * Γ(X_t) / γA
ωB(X_t) = (1 - X_t) * Γ(X_t) / γB;
$$\Delta(X_t) = \sum_i \omega^i(X_t)\frac{\gamma^i + 1}{\gamma^i}$$
Δ(X_t) = ωA(X_t) * (γA + 1) / γA + ωB(X_t) * (γB + 1) / γB;
The primitive constants of the model are:
const vA = 0.01; const vB = 1 - vA;
const ρ = 0.001
const δ1 = 0.0525; const δ2 = 0.0611
const ω = 1 - 0.08
const γA = 1.5; const ψA = 0.70; const αA = 1 - 1/ψA
const γB = 10.; const ψB = 0.05; const αB = 1 - 1/ψB
const μ_Y = 0.02; const σ_Y = 0.041
const μ = 0.02; const σ = 0.041
const π = 0.02;
Applying the normalization of the integral on the right-hand side in equation (4) in the paper to one yields:
const B1 = 30.72 / (π / (π + δ1) * 30.72 + π / (π + δ2) * -30.29)
const B2 = -30.29 / (π / (π + δ1) * 30.72 + π / (π + δ2) * -30.29);
The auxiliary constants are defined by:
$M_1^i = - 1 - \frac{\Xi_1^i}{\gamma^i}$, $M_2^i = \frac{\gamma^i - 1}{\gamma^i}$, $\Xi_1^i = -\frac{\alpha^i+\gamma^i-1}{\alpha^i}$, $\Xi_2^i = \frac{\alpha^i}{(1-\alpha^i)(1-\gamma^i)}$, $\Xi_3^i = -\frac{\rho + \pi}{\alpha^i}(1-\gamma^i)$, $\Xi_4^i = - \frac{\alpha^i + \gamma^i - 1}{(1-\alpha^i)(1-\gamma^i)}$
const ΞA_1 = - (αA + γA - 1) / αA; const ΞB_1 = - (αB + γB - 1) / αB
const ΞA_2 = αA / ((1-αA) * (1-γA)); const ΞB_2 = αB / ((1-αB) * (1-γB))
const ΞA_3 = - (ρ + π) / αA * (1 - γA); const ΞB_3 = - (ρ + π) / αB * (1 - γB)
const ΞA_4 = - (αA + γA - 1) / ((1-αA) * (1-γA)); const ΞB_4 = - (αB + γB - 1) / ((1-αB) * (1-γB))
const MA_1 = - 1 - ΞA_1 / γA; const MB_1 = - 1 - ΞB_1 / γB
const MA_2 = (γA - 1) / γA; const MB_2 = (γB - 1) / γB;
The following function carries out finite differencing. It takes the evaluation of a solution function on a grid as input and calculates the first derivative via upwinding and the second derivative via central differences. Hence, f, fx, fxx, and wind are vectors with length N.
function diff!(f, fx, fxx, wind, N)
for i in 1:N
#fx via upwinding
if wind[i] >= 0
fx[i] = (f[min(i + 1, N)] - f[i]) * (N - 1)
else
fx[i] = (f[i] - f[max(i - 1, 1)]) * (N - 1)
end
#fxx not via upwinding
fxx[i] = (f[min(i + 1, N)] * (N - 1)^2 + f[max(i - 1, 1)] * (N - 1)^2 - 2 * f[i] * (N - 1)^2)
end
end
Now, the right-hand side of equations (1) for $i\in\{A,B\}$ and equation (2) for $j\in\{1,2\}$ are implemented. The following function takes a stacked candidate solution of the four solution functions evaluated on a grid as input and returns the residual of those four equations.
It carries out the following steps:
function F(Y)
#unstack the input into the four solution functions
N = round(Int64, length(Y)/4)
x = collect(linspace(0, 1, N))
ϕ1 = Y[1:N]
ϕ2 = Y[N+1:2*N]
pA = Y[2*N+1:3*N]
pB = Y[3*N+1:4*N]
#setup derivatives
ϕ1x = similar(ϕ1)
ϕ1xx = similar(ϕ1)
ϕ2x = similar(ϕ2)
ϕ2xx = similar(ϕ2)
pAx = similar(pA)
pAxx = similar(pA)
pBx = similar(pB)
pBxx = similar(pB)
#finite differences without upwinding
diff!(ϕ1, ϕ1x, ϕ1xx, zeros(200), N)
diff!(ϕ2, ϕ2x, ϕ2xx, zeros(200), N)
diff!(pA, pAx, pAxx, zeros(200), N)
diff!(pB, pBx, pBxx, zeros(200), N)
#update the wind
wind = μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2)
diff!(ϕ1, ϕ1x, ϕ1xx, wind, N)
diff!(ϕ2, ϕ2x, ϕ2xx, wind, N)
diff!(pA, pAx, pAxx, wind, N)
diff!(pB, pBx, pBxx, wind, N)
#the ODEs, stacked into one single vector
vcat(
1/2 * ϕ1xx .* σ_X.(x, pA, pAx, pB, pBx).^2 +
ϕ1x .* (μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + σ_X.(x, pA, pAx, pB, pBx) .* (σ_Y - κ.(x, pA, pAx, pB, pBx))) +
ϕ1 .* (μ_Y - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - π - δ1 - σ_Y .* κ.(x, pA, pAx, pB, pBx)) +
B1 * ω,
1/2 * ϕ2xx .* σ_X.(x, pA, pAx, pB, pBx).^2 +
ϕ2x .* (μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + σ_X.(x, pA, pAx, pB, pBx) .* (σ_Y - κ.(x, pA, pAx, pB, pBx))) +
ϕ2 .* (μ_Y - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - π - δ2 - σ_Y .* κ.(x, pA, pAx, pB, pBx)) +
B2 * ω,
pA .* (1 ./ pA +
ψA * (r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - ρ) + κ.(x, pA, pAx, pB, pBx).^2 * (1 + ψA) / (2 * γA) + (1 - ψA * γA) / (γA * (ψA - 1)) * κ.(x, pA, pAx, pB, pBx) .* (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) -
(1 - γA * ψA) / (2 * γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)).^2 - π + pAx ./ pA .* μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + 0.5 * pAxx ./ pA .* σ_X.(x, pA, pAx, pB, pBx).^2 +
(κ.(x, pA, pAx, pB, pBx) / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx))) .* (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - κ.(x, pA, pAx, pB, pBx) .* ((pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) + (κ.(x, pA, pAx, pB, pBx) / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx))))),
pB .* (1 ./ pB +
ψB * (r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - ρ) + κ.(x, pA, pAx, pB, pBx).^2 * (1 + ψB) / (2 * γB) + (1 - ψB * γB) / (γB * (ψB - 1)) * κ.(x, pA, pAx, pB, pBx) .* (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) -
(1 - γB * ψB) / (2 * γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)).^2 - π + pBx ./ pB .* μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + 0.5 * pBxx ./ pB .* σ_X.(x, pA, pAx, pB, pBx).^2 +
(κ.(x, pA, pAx, pB, pBx) / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx))) .* (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - κ.(x, pA, pAx, pB, pBx) .* ((pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) + (κ.(x, pA, pAx, pB, pBx) / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)))))
)
end
One approach is to solve $F(Y) = 0$ by means of guessing an initial value for $Y$ and finding $Y_{t+1}$ from:
$$0=F(Y_t) + J_F(Y_t)\Big(Y_{t+1} - Y_t\Big)$$
This usually needs a good starting value, though. Gârleanu and Panageas solve the ODE with this approach, but use a good starting value.
The following function starts with an initial guess for $Y$ and iterates on the above equation, by calculating the Jacobian via automatic differentiation:
using ForwardDiff
function ExplicitTimeStepping(;initial_Y = ones(800), verbose = true)
Yold = initial_Y
Ynew = zeros(Yold)
distance = 0.1
iteration = 0
while (distance > 1e-10) && (iteration < 20)
Ynew .= -ForwardDiff.jacobian(F, Yold)^-1 * F(Yold) + Yold
distance = maximum(abs.(Yold .- Ynew))
Yold = copy(Ynew)
iteration += 1
if verbose && (iteration % 1 == 0)
println("Iteration $iteration with distance $distance")
end
end
return Ynew
end
ExplicitTimeStepping();
This is highly unstable. The starting value is not good enough!
Another approach is to solve $F(Y) = 0$ by means of guessing an initial value for $Y$ and finding $Y_{t+1}$ from:
$$0=F(Y_{t+1}) + \frac{1}{\Delta}\Big(Y_{t+1} - Y_t\Big)$$
where $\Delta$ is the step size.
Matthieu Gomez suggests solving this equation via a Newton-Ralphson method in his package EconPDEs, see his notes.
The proposed algorithm has two nested iterations:
The outer iteration over $t$ only concerns the $Y_t$. The outer iteration:
The inner iteration takes $Y_t$ and the step size as given. The above equation is solved exactly for $Y_t+1$ by a Newton-Ralphson method. That is, it iterates over $i$ and begin with $Y_{t+1}^0 = Y_t$ as starting value:
$$ 0 =F(Y_{t+1}^{i})+\frac{1}{\Delta}\Big(Y_{t+1}^{i}-Y_{t}\Big)+ \Big(J_{F}(Y_{t+1}^{i})-\frac{1}{\Delta}\Big)\Big(Y_{t+1}^{i+1}-Y_{t+1}^{i}\Big) $$
This is the first-order Taylor approximation of the above equation expanded around $Y^i_{t}$ and it can be easily solved for $Y^i_{t+1}$. Mind, the Jacobian is a banded matrix, which could be exploited for higher computational performance.
The following function encloses $Y_t$ and the step size. That is, it returns an anonymous function that knows the values of those two objects. The returned function takes a candidate solution of the inner iteration Yprime and an residual vector as input and calculates the residual:
function residual_wrapper(Y, stepsize = 1.)
return (residual, Yprime) -> residual .= F(Yprime) - 1/stepsize * (Yprime - Y)
end
The following function conducts the two nested iteration. The inner iteration is simply a call of the nlsolve function with method = :newton.
using NLsolve
function ImpicitTimeStepping(;initial_Y = ones(800), initial_stepsize = 1., verbose = true)
Yold = initial_Y
Ynew = zeros(Yold)
stepsize = initial_stepsize
distance = 0.1
distanceold = 0.1
iteration = 0
while (distance > 1e-10) && (iteration < 200) && (stepsize >= 1e-12)
iteration += 1
result = nlsolve(residual_wrapper(Yold, stepsize), Yold, iterations = 25,
autodiff=:forward, method = :newton, ftol = 1e-9)
Ynew .= result.zero
if any(isnan.(Ynew))
println("Iteration $iteration bad solution with step size $stepsize")
stepsize = stepsize / 10
continue
end
if !result.f_converged
println("Iteration $iteration no convergence with step size $stepsize")
stepsize = stepsize / 10
continue
end
distance, distanceold = maximum(abs.(F(Ynew))), distance
Yold = copy(Ynew)
if (distance <= distanceold)
stepsize = stepsize * 10
end
if verbose && (iteration % 1 == 0)
println("Iteration $iteration with distance $distance")
end
end
return Ynew
end
Despite using a naive starting value of 800 ones, the algorithm performs very well:
@time Y = ImpicitTimeStepping();
In comparison to the explicit time-stepping method, starting with a candidate that is less than one percent off from the solution maximally (according to the discretization) fails to converge sometimes:
# scatter plot↔
# plot solution functions↔
The following figure replicates Figure 1 from the paper. The volatility of stock prices is different, which yields a different equity premium. This may be a result of the error in equation (A.16) in the paper.
# plot other functions↔
# plot the law of motion for X↔
Now, the process for $X_t$ is simulated in order to find its stationary distribution:
# plot long path of X↔
Table 1 in the paper shows the assets pricing relevant objects, which are calculated here again:
# draw table↔
It appears that the model implied parameters are lower than reported in the paper.